\section{多元线性回归习题}

\textbf{3} 证明$\hat{\sigma}^{2}=\frac{1}{n-p-1}SSE$是误差项方差$\sigma^{2}$的无偏估计。

\begin{proof}[\bf 解]

\begin{equation}
    \label{eq_sigma2}   
    \hat{\sigma}^{2} = \frac{1}{n-p-1} \sum_{i=1}^{n} e_{i}^{2}  
\end{equation}

对式\eqref{eq_sigma2}两边同时取期望，得到：

\begin{equation}
    E(\hat{\sigma}^{2}) = E[\frac{1}{n-p-1}\sum_{i=1}^{n}e_{i}^{2}]
\end{equation}
由于$\epsilon_{i}$之间相互独立,且$E(\epsilon_{i})=0$，故
\begin{eqnarray}
     E(\hat{\sigma}^{2}) &=& \frac{1}{n-p-1}E[\sum_{i=1}^{n}e_{i}^{2}] \nonumber \\
     &=& \frac{1}{n-p-1} \sum_{i=1}^{n} Var(e_{i})  \label{Var_e2}
\end{eqnarray}
故核心问题为求残差$e_{i}$的方差。根据残差的定义有
\begin{eqnarray}
    \boldsymbol{e} &=& \boldsymbol{y} - \hat{\boldsymbol{y}} \nonumber \\
    & = &  \boldsymbol{y} - X\hat{\beta} \nonumber \\
    & = &  \boldsymbol{y} - X(X^{\top}X)^{-1}X^{\top}\boldsymbol{y} \nonumber \\
    & = &  (I - H)\boldsymbol{y}
\end{eqnarray}
其中$I$为单位矩阵，$H$为帽子矩阵$H = X(X^{\top}X)^{-1}X^{\top}$。
\begin{eqnarray}
    Var(\boldsymbol{e}) &=& Cov(\boldsymbol{e},\boldsymbol{e}) \nonumber \\
    & = & Cov((I-H)\boldsymbol{y},(I-H)\boldsymbol{y}) \nonumber \\
    & = & (I-H)Cov(\boldsymbol{y},\boldsymbol{y})(I-H)^{\top} \nonumber \\
    & = & \sigma^{2}(I-H)I_{n}(I-H) \nonumber \\
    & = & \sigma^{2}(I-H) \nonumber
\end{eqnarray}
易得
\begin{equation}
    Var(e_{i}) = \sigma^{2}(1 - h_{ii}), i=1,2,\cdots,n \label{vare2}
\end{equation}
将\eqref{Var_e2}代入\eqref{vare2}，可得
\begin{equation}
    E(\hat{\sigma^{2}}) = \frac{1}{n-p-1} \sigma^{2} \sum_{i}^{n} (1-h_{ii}) = \frac{1}{n-p-1} \sigma^{2} (n - \sum_{i}^{n}h_{ii}) \label{sigma2_3}
\end{equation}
其中$\sum_{i}^{n}h_{ii} = tr(H)$,故求解帽子矩阵的迹
\begin{eqnarray}
    tr(H) &=& tr(X(X^{\top}X)^{-1}X^{\top}) \nonumber \\
          &=& tr((X^{\top}X)^{-1}XX^{\top}) \nonumber \\
          &=& tr(I_{p+1}) = p+1 \nonumber
\end{eqnarray}
即
\begin{equation}
    \sum_{i=1}^{n}h_{ii} = p+1
\end{equation}
代入\eqref{sigma2_3}可得
\begin{equation}
    E(\hat{\sigma}^{2}) = \frac{1}{n-p-1} \sigma^{2} (n - p - 1) = \sigma^{2} 
\end{equation}
\end{proof}

\textbf{4}一个回归方程的复相关系数$R=0.99$,样本决定系数$R^{2} = 0.9801$，我们能断定回归方程很理想吗？
\begin{proof}[\textbf{答}]
不能。$R^{2}$的大小与自变量的数目及样本量有关。当样本量$n$与自变量的个数接近时，$R^{2}$容易接近$1$，这其中隐含着一些虚假的成分。并且通过增加自变量的个数会使得$R^{2}$只增加而不会减少，但是过多的无关变量会使得估计的方差变大，故$R^{2}$仅能作为模型拟合程度的一个参考指标。
\end{proof}

\textbf{7}验证式$\hat{\beta}^{*}_{j} = \frac{\sqrt{L_{jj}}}{\sqrt{L_{yy}}}\hat{\beta_{j}},j=1,2,\cdots,p$
\begin{proof}
原始经验回归方程为：
\begin{equation}
    \label{emp_eq}
    \hat{y}_{i} = \hat{\beta}_{0} + \hat{\beta}_{1}x_{i 1} + \hat{\beta}_{2}x_{i 2} + \cdots \hat{\beta}_{p}x_{i p}.
\end{equation}
由经验回归方程过样本均值可得，
\begin{equation}
    \label{e_eq}
    \bar{y} = \hat{\beta}_{0} + \hat{\beta}_{1}\bar{x}_{1} + \hat{\beta}_{2}\bar{x}_{2} + \cdots \hat{\beta}_{p}\bar{x}_{p}.
\end{equation}

使用式\eqref{emp_eq}减去式\eqref{e_eq}得到

\begin{equation}
    \label{center_eq}
    \hat{y}_{i} - \bar{y} = \hat{\beta}_{1}(x_{i 1} - \bar{x}_{1}) + \hat{\beta}_{2}(x_{i 2} - \bar{x}_{2}) + \cdots \hat{\beta}_{p}(x_{i p} - \bar{x}_{p}).
\end{equation}
对式\eqref{center_eq}两边同时除以$\sqrt{L_{y y}}$，并对其变形整理得到

\begin{equation}
    \frac{\hat{y}_{i} - \bar{y}}{\sqrt{L_{y y}}} = \frac{\hat{\beta}_{1}(x_{i 1} - \bar{x}_{1})}{\sqrt{L_{1 1}}}\times \frac{\sqrt{L_{1 1}}}{\sqrt{L_{y y}}} + \frac{\hat{\beta}_{2}(x_{i 2} - \bar{x}_{2})}{{\sqrt{L_{2 2}}}}\times \frac{\sqrt{L_{2 2}}}{\sqrt{L_{y y}}} + \cdots \frac{\hat{\beta}_{p}(x_{i p} - \bar{x}_{p})}{{\sqrt{L_{3 3}}}} \times \frac{\sqrt{L_{p p}}}{\sqrt{L_{y y}}}.
\end{equation}
将$\hat{y}^{*}_{i} = \frac{\hat{y}_{i} - \bar{y}}{\sqrt{L_{y y}}} , x^{*}_{i j} = \frac{x_{i j}-\bar{x_{i}}}{\sqrt{L_{j j}}}$代入上式可得
\begin{equation}
    \hat{y}^{*}_{i} = \frac{\sqrt{L_{11}}}{\sqrt{L_{y y}}}\hat{\beta}_{1}x_{i}^{*} + \frac{\sqrt{L_{2 2}}}{\sqrt{L_{y y}}}\hat{\beta}_{2}x_{i}^{*} + \cdots + \frac{\sqrt{L_{p p}}}{\sqrt{L_{y y}}}\hat{\beta}_{p}x_{i}^{*}.  
\end{equation}
即
\begin{equation}
    \hat{\beta}^{*}_{j} = \frac{\sqrt{L_{jj}}}{\sqrt{L_{yy}}}\hat{\beta_{j}},j=1,2,\cdots,p. \nonumber
\end{equation}
\end{proof}

\textbf{8}利式（3.60）证明式（3.61）成立，即
\begin{equation}
    r_{12,3} = \frac{r_{12}-r_{13}r_{23}}{(\sqrt{1-r_{13}^{2})(\sqrt{1-r_{23}^{2}})}}. \nonumber
\end{equation}

\begin{proof}
样本的相关系数矩阵为:
\begin{equation}
r = \left|\begin{matrix}
1      & r_{12} & r_{13}\\
r_{12} & 1      & r_{23}\\
r_{13} & r_{23} & 1
\end{matrix}\right|
\end{equation} 
再根据式(3.60)
\begin{equation}
    \label{eq_23}
    r_{12,3} = \frac{-\Delta_{12}}{\sqrt{\Delta_{11}\Delta_{22}}}.
\end{equation}
其中
\begin{equation}
-\Delta_{12} = -(-1)^{1+2} 
\left|\begin{matrix}
r_{12} &  r_{23}\\
r_{13} & 1
\end{matrix}\right| = r_{12} - r_{13}r_{23} \nonumber
\end{equation}
\begin{equation}
\Delta_{11} = (-1)^{1+1} 
 \left|\begin{matrix}
1      &  r_{23}\\
r_{23} & 1
\end{matrix}\right| = 1-r_{23}^{2} \nonumber
\end{equation}

\begin{equation}
\Delta_{22} = (-1)^{2+2} 
\left|\begin{matrix}
1      &  r_{13}\\
r_{13} & 1
\end{matrix}\right| = 1-r_{13}^{2} \nonumber
\end{equation}
带入\eqref{eq_23}可得
\begin{equation}
r_{12,3} = \frac{r_{12}-r_{13}r_{23}}{(\sqrt{1-r_{13}^{2})(\sqrt{1-r_{23}^{2}})}}. \nonumber \nonumber\end{equation}
\end{proof}

\textbf{10} 验证决定系数$R^{2}$与$F$值之间的关系式
\begin{equation}
    R^{2} = \frac{F}{F + (n-p-1)/p} \nonumber
\end{equation}
\begin{proof}

根据$F$统计量的表达式：
\begin{equation}
    F = \frac{SSR/p}{SSE/(n-p-1)} \nonumber
\end{equation}
可得
\begin{equation}
    \label{ch_2_eq_25}
    SSR = \frac{pF\times SSE}{n-p-1}
\end{equation}
再由决定系数表达式
\begin{equation}
    \label{ch_2_eq26}
    R^{2} = \frac{SSR}{SST} = \frac{SSR}{SSE+SSR}
\end{equation}
将式\eqref{ch_2_eq_25}代入式\eqref{ch_2_eq26}可得：
\begin{eqnarray}
    R^{2} &=& \frac{\frac{pF\times SSE}{n-p-1}}{SSE+\frac{pF\times SSE}{n-p-1}} \nonumber\\
    & = & \frac{\frac{pF}{n-p-1}}{1+\frac{pF}{n-p-1}}  \nonumber \\
    & = & \frac{pF}{n-p-1 + pF} \nonumber \\
    & = & \frac{F}{F + (n-p-1)/p}  \nonumber
\end{eqnarray}
    
\end{proof}

\textbf{11}研究货运总量$y$(万吨)与工业总产值$x_{1}$(亿元)、农业总产值$x_{2}$、农业总产值$x_{2}$(亿元)之间的关系。数据如表\ref{tab_2_11}所示

\begin{table}[H]
\caption{回归数据}
\label{tab_2_11}
\centering
\begin{tabular}{@{}ccccc@{}}
\toprule
\toprule
编号 & 货运总量y(万吨) & 工业总产值$x_{1}$(亿元) & 农业总产值$x_{2}$(亿元) & 居民非商品支出$x_{3}$(亿元) \\ \midrule
1  & 160       & 70               & 35               & 1.0                \\
2  & 260       & 75               & 40               & 2.4                \\
3  & 210       & 65               & 40               & 2.0                \\
4  & 265       & 74               & 42               & 3.0                \\
5  & 240       & 72               & 38               & 1.2                \\
6  & 220       & 68               & 45               & 1.5                \\
7  & 275       & 78               & 42               & 4.0                \\
8  & 160       & 66               & 36               & 2.0                \\
9  & 275       & 70               & 44               & 3.2                \\
10 & 250       & 65               & 42               & 3.0                \\ \bottomrule
\end{tabular}
\end{table}

(1)计算出$y,x_{1},x_{2},x_{3}$的相关系数矩阵。

计算相关系数的R代码如下：
\begin{lstlisting}[language=R]
    > x1 <- c(70, 75, 65, 74, 72, 68, 78, 66, 70, 65)
    > x2 <- c(35, 40, 40, 42, 38, 45, 42, 36, 44, 42)
    > x3 <- c(1.0, 2.4, 2.0, 3.0, 1.2, 1.5, 4.0, 2.0, 3.2, 3.0)
    > y  <- c(160, 260, 210, 265, 240, 220, 275, 160, 275, 250)
    > dataset <- data.frame(y=y,x1=x1,x2=x2,x3=x3)
    > cor(dataset)

           y        x1        x2        x3
    y  1.0000000 0.5556527 0.7306199 0.7235354
    x1 0.5556527 1.0000000 0.1129513 0.3983870
    x2 0.7306199 0.1129513 1.0000000 0.5474739
    x3 0.7235354 0.3983870 0.5474739 1.0000000
\end{lstlisting}

(2)求$y$关于$x_{1},x_{2},x_{3}$的三元线性回归方程。
model = lm(y~x1+x2+x3)
model
\begin{lstlisting}[language=R]
    > model = lm(y~x1+x2+x3)
    > model
    Call:
    lm(formula = y ~ x1 + x2 + x3)
    
    Coefficients:
    (Intercept)           x1           x2           x3  
      -348.280          3.754         7.101       12.447  
\end{lstlisting}
$y$关于$x_{1},x_{2},x_{3}$的三元线性回归方程为：
\begin{equation}
    y = 3.754x_{1} + 7.101x_{2} + 12.447x_{3} -348.280. \nonumber
\end{equation}

(3) 对所得方程做拟合优度检验。

根据代码结果，模型的$R^{2}=0.8055$拟合优度较好。
\begin{lstlisting}
    > summary(model)
    Call:
    lm(formula = y ~ x1 + x2 + x3)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -25.198 -17.035   2.627  11.677  33.225 
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
    (Intercept) -348.280    176.459  -1.974   0.0959 .
    x1             3.754      1.933   1.942   0.1002  
    x2             7.101      2.880   2.465   0.0488 *
    x3            12.447     10.569   1.178   0.2835  
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 23.44 on 6 degrees of freedom
    Multiple R-squared:  0.8055,	Adjusted R-squared:  0.7083 
    F-statistic: 8.283 on 3 and 6 DF,  p-value: 0.01487
\end{lstlisting}



(4) 对回归方程做显著性检验。

根据(3)中的结果，$F$统计量的观测值为8.283，对应的概率$P$值为0.01487在显著性水平$\alpha =5\%$的条件下回归方程显著，即回归系数不全为0。

(5) 对每一个回归系数做显著性检验。

根据(3)中的结果，$x_{1}，x_{2}$的系数可认为在10\%水平上显著，而$x_{3}$的P值为0.2835不显著。

(6) 如果有的回归系数没通过显著性检验，将其剔除，重新简历回归方程，再做方程的显著性检验和回归系数的显著性检验。

剔除$x_{3}$
\begin{lstlisting}
    > model_2 <- lm(y~x1+x2)
    > summary(model_2)
    Call:
    lm(formula = y ~ x1 + x2)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -42.012 -10.656   4.358  11.984  28.927 
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)   
    (Intercept) -459.624    153.058  -3.003  0.01986 * 
    x1             4.676      1.816   2.575  0.03676 * 
    x2             8.971      2.468   3.634  0.00835 **
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 24.08 on 7 degrees of freedom
    Multiple R-squared:  0.7605,	Adjusted R-squared:  0.6921 
    F-statistic: 11.12 on 2 and 7 DF,  p-value: 0.006718
\end{lstlisting}

根据代码运行的结果，$F$统计量的观测值为11.12对应的概率$P$值为0.006718，在显著性水平$\alpha=0.01$的水平上显著，即回归方程的系数不全为0。对于每个回归系数，都在5\%水平通过显著性检验。

(7) 求出每一个回归系数的95\%置信区间。

\begin{lstlisting}
    > confint(model,level = 0.95)
                    2.5 %    97.5 %
    (Intercept) -780.0603287 83.499990
    x1            -0.9766149  8.484688
    x2             0.0529179 14.148507
    x3           -13.4147488 38.309689
\end{lstlisting}

故$\beta_{1},\beta_{2},\beta_{3}$的置信区间分别为:[-0.9766149,8.484688],[0.0529179,14.148507],[-13.4147488, 38.309689]

(8) 求出标准化回归方程

\begin{lstlisting}
    > z_x1 = scale(x1)
    > z_x2 = scale(x2)
    > z_x3 = scale(x3)
    > z_y = scale(y)
    > model = lm(z_y~z_x1+z_x2+z_x3)
    > model

    Call:
    lm(formula = z_y ~ z_x1 + z_x2 + z_x3)
    
    Coefficients:
    (Intercept)         z_x1         z_x2         z_x3  
     -6.063e-16    3.848e-01    5.355e-01    2.771e-01 
\end{lstlisting}

标准化回归方程为：
\begin{equation}
    y^{*} = 0.385x_{1}^{*} + 0.536x_{2}^{*} + 0.277x_{3}^{*}
\end{equation}
(9) 求当$x_{01}=75,x_{02}=42,x_{03}=3.1$时的$\hat{y}_{0}$，给定置信水平为95\%，计算精确置信区间，并且手工计算近似预测区间。


\begin{lstlisting}
    > x0 = data.frame('x1'=75,'x2'=42,'x3'=3.1)
    > predict(model, newdata = x0,interval ="prediction")
           fit      lwr      upr
    1 270.0897 206.059 334.1204
\end{lstlisting}

$\hat{y}_{0}=270.0897$ 

精确预测区间为：[206.059,334.1204]

根据(3)，$\hat{\sigma} = 23.44$，故近似预测区间为：[223.2097, 316.9697]

(10) 结合回归方程对问题做一些基本分析。

根据回归方程的结果方程整体显著，但在其中$x_{3}$的系数不显著，说明$x_{3}$对y没有显著的作用，整体拟合程度较好$R^{2}=0.8055$。